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ABSTRACT 

We study the short- and long-term effects of an intermediate mass black hole (IMBH) on the orbits of stars 
bound to the supermassive black hole (SMBH) at the center of the Milky Way. A regularized A^-body code 
including post-Newtonian terms is used to carry out direct integrations of 19 stars in the S-star cluster for 
lOMyr. The mass of the IMBH is assigned one of four values from 400 M© to 4000 M©, and its initial semi- 
major axis with respect to the SMBH is varied from 0.3-30mpc, bracketing the radii at which inspiral of the 
IMBH is expected to stall. We consider two values for the eccentricity of the IMBH/SMBH binary, e = (0,0.7), 
and 12 values for the orientation of the binary's plane. Changes at the level of ~ 1% in the orbital elements of 
the S-stars could occur in just a few years if the IMBH is sufficiently massive. On time scales of 1 Myr or longer, 
the IMBH efficiently randomizes the eccentricities and orbital inclinations of the S-stars. Kozai oscillations are 
observed when the IMBH lies well outside the orbits of the stars. Perturbations from the IMBH can eject 
stars from the cluster, producing hypervelocity stars, and can also scatter stars into the SMBH; stars with high 
initial eccentricities are most likely to be affected in both cases. The distribution of S-star orbital elements 
is significantly altered from its currently-observed form by IMBHs with masses greater than ~ 10^ Mq if the 
IMBH/SMBH semi-major axis lies in the range 3- lOmpc. We use these results to further constrain the allowed 
parameters of an IMBH/SMBH binary at the Galactic center. 
Subject headings: Galaxyxenter - stellar dynamics 



1. INTRODUCTION 

Near-infrared imaging and spectroscopy of the Galactic 
center indicate the presence of over 100 young massive stars 
in the inner parsec of the Milky Way. These belong to two dis- 
tinct populations: (1) Outside the central arcsecond, a group 
of roughly 40 stars, mostly O and Wolf-Rayet stars with an 
estimated age of 6 Myr, move on approxim ately circular or- 
bits in a thin and coherent stellar disk (Pau mard et al.ll2006h 
(2) Inside the central arcsecond, a group of 20 B-type stars 
(the S-stars) move on eccentric and randomly oriented or- 
bits around the supermassive black hole (SMBH) dGhez et al.l 
|2003; Eisenhauer et al. 2005; Gillessen et al. 2009^ 

The existence of young massive stars so close to the SMBH 
is puzzling since tidal stresses from the SMBH are likely to 
suppress star formation in this region. If the S-stars formed 
at larger distances from Sgr A* and were then transported in- 
ward, the migration would need to occur on a timescale of 
only a few Myr, given the apparent age of the stars. This con- 
stitutes the so-called paradox of youth. 

Several scenarios have been proposed for the origin of the 
young stars in the Galactic center, including (i) rejuvenation 
of older stars due to physical col lisions and/or tid al stripping 
in the vicinity of the SMBH (?Ghez et al. 2003), (ii) tidal 
disruption of binary stars on low angular momentum orbits 
(iGouId & O uillen 2003; ?), (iii) in-situ formation from a frag- 
men ting gas disk forme d from an infalling molecular cloud 
(?Bo nnell & Ricell2008l) . (iv) formation in a star cluster out- 
side the central parsec and transportation by the inspiralling 
cluster (Gerhard 2001). 

The spectra of the observed S-stars seem to indicate that 
they are ordinary B-type main-sequ ence stars, with no ap- 
parent evide nce for an exotic history (lEisenhauer et al.lf2005t 
lFigeiil2008l) . This would seem to exclude the rejuvenation 
scenario. 

Electronic address: [alessiag,merritt@astro.rit.edu] 



The binary disruption scenario postulates that the S-stars 
were once in binary systems on plunging orbits that were 
tidally disrupted by the SMBH, ejecting one star and leaving 
the other in an extreme ly eccentric orbit around the SMBH 
dGould & Ouillerj2003h . This model requires an ad hoc reser- 
voir of new binaries at large radii, a mechanism to scatter 
the stars onto p lunging orbits (e.g. "massive perturbers"; 
iPerets etaP ('2007)), and an efficient mechanism to thermalize 
the eccentricities within the lifetime of the stars (e.g. "reso- 
nant r elaxation" from stellar mass black holes; iPerets etaP 
(120091) ). Modified binary disruption scenarios in which the 
S-stars were born in a stellar disk have also been proposed 
(??). 

The in-situ formation model assumes that an infalling 
molecular cloud interacts with the SMBH and forms a disk 
which then fragments and forms stars. This model success- 
fuUy reproduces the observed properties of the disk stars 
dBonnell & Ricd 12008 ), including mild eccentricities and a 
top-heavy mass function. It however cannot account for the 
large eccentricities and random inclinations of the orbits of 
th e S-stars. 

iMerritt et al.l ([2009) showed that these dynamical proper- 
ties can be naturally explained by interaction of the stars with 
an intermediate mass black hole (IMBH) in the context of 
the inspiralling star cluster scenario. This model postulates 
that the S-stars were born in a star cluster at large enough 
distance from the SMBH that a giant molecular cloud could 
collapse and fragment. The new cluster then migrated in- 
ward due to dynamical friction, depositing stars along the 
way as the tidal force from the SMBH stripped them from 
the outer regions of the cluster (Gerhard 200 JJ. If the clus- 
ter has by then formed an IMBH at its center, e.g. via 
the ru naway merging scenario ( Portegies Zwart & McMillaiil 
l2002h . the cluster can remain bound and spiral inward 
for much longer before releasing all its component stars 
dHansen & Milosavl)evicll2003.) . reaching distances of a few 
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tens of milliparsecs before evolutio n of the IMBH/SMBH bi- 
nary stalls ( Baumgardt et al. 2006; Matsubavashi et ai]|2007l : 
iLQckmann & Baumgardti2008 ). Interactions with the same 
IMBH that brought the stars close to the SMBH can then 
randomize their orbital inclinations and thermalize th eir ec- 
centri city distribution on timescales of 1 Myr (Merri tt et alJ 
l2009h . The IMBH can also capture stars ejected from the clus- 
ter into a resonance and bring them close to the SMBH (?). 

Whether or not the inspiralling star cluster model is cor- 
rect, this work raises intriguing questions about the observ- 
able consequences of an IMBH that is currently imbedded in 
the S-star cluster. Are perturbations from an IMBH detectable 
on short (year-long) time scales, and if so, which of the S- 
stars are most likely to be affected? How stable is the cluster 
on time scales of order stellar lifetimes? Would an apprecia- 
ble fraction of the stars be ejected from the Galactic center, or 
scattered into the SMBH, by the IMBH? What constraints can 
be placed on the physical parameters (mass, orbital elements) 
of a putative IMBH/SMBH binary given the need to conserve 
the observed distribution of orbital elements of the S-stars? 

To address these questions, we have carried out a large 
set of high-accuracy A^-body simulations, taking as our ini- 
tial conditions the current positions and velocities of the ~ 20 
S-stars with well determined orbits, plus an additional mas- 
sive particle representing the IMBH. In §2 we describe the 
initial conditions and the A^-body methods. §3 discusses the 
short (few years) time scale effects on the observed stellar or- 
bits. The longer term (Myr) behavior of the S-star cluster is 
described in §4. §5 summarizes our results and draws conclu- 
sions about the allowed parameter space for an IMBH/SMBH 
binary at the Galactic center 

2. INITIAL MODELS AND NUMERICAL METHODS 

Our A^-body mo dels include a SMBH, a n IMBH and a clus- 
ter of 19 S-stars. iGillessen et al.l ( 12009 *) provide Keplerian 
elements for 28 stars, of which one (Sill) appears to be un- 
bound as a result of its large radial velocity, six (S66, S67, 
S83, S87, S96, S97) likely belong to the clockwise stellar disk 
dGenzel et al.ll2003HPaumard et al.ll2 006) and one (S71) has a 
very long orbital period. The remaining stars have well de- 
fined or at least reliable orbits and exhibit essentially random 
orbital orientations. We determined positions and velocities 
at year 2008 AD for the 19 remaining stars from the classi- 
cal elements' given in Table[ T| The assu med SMBH mass is 
MsMBH = 4.0 X IO^Mq (iGillessen et al. 200 9J. The masses 
of the S-stars were set to IOMq (e.g. Eisenh auer et alj|2005[) 
except for s tar S-2 for which a value of 20 M0 was adopted 
dMartins et al. 2008). 

We considered four values of the IMBH mass, 
MiMBH = (400, 1000, 2000, 4000) Mq, or ^ = MsmbhMmbh = 
(1,2.5,5, 10) X lO"''; five values of the initial semi-major 
axis, a = (0.3, 1,3, 10,30)mpc; two values of the eccentricity, 
e = (0,0.7); and twelve orientations of the IMBH/SMBH 
orbital plane, for a total of 480 initial models. The parameters 
of the simulations are summarized in TablelH 

' The orbit al elements are defined as in lEisenh auer et all )200 5D and 
IPaumard et al] ( 2006): a = semi-major axis; P = orbital period; e = eccen- 
tricity; / = inclination of the ascending node, i.e. angle between the orbital 
plane and the plane tangential to the celestial sphere (0° < < 90° for direct 
counterclockwise projected rotation, 90° < r < 1 80° for retrograde clockwise 
projected rotation); f! = position angle of the ascending node; lu = longitude 
of periastron, i.e. the angle between the radius vector of the ascending node 
and that of the periastron counted from the node in the direction of the orbital 
motion. 



TABLE 1 
S-STAR ORBITAL ELEMENTS 



star 


^ orb 


£1 


f> 


i 


SI 












(dee) 






SI 


128.0 


19.462 


0.483 


120.1 


341.3 


116.9 


S2 


15.8 


4.829 


0.880 


134.5 


225.8 


63.8 


S4 


58.2 


11.523 


0.407 


77.8 


257.9 


316.9 


S5 


210.0 


27.072 


0.672 


1 10.2 


147.0 


270.6 


S6 


119.0 


18.601 


0.898 


86.5 


84.0 


126.2 


S8 


95.6 


16.038 


0.816 


74.6 


315.7 


344.4 


S9 


62.3 


12.054 


0.867 


80.4 


145.7 


229.7 


S12 


60.1 


11.763 


0.899 


32.8 


236.9 


311.2 


S13 


56.3 


11.263 


0.471 


27.8 


72.0 


247.5 


S14 


45.2 


9.727 


0.962 


99.4 


227.9 


339.0 


S17 


62.0 


12.019 


0.370 


96.3 


187.2 


318.7 


S18 


47.9 


10.123 


0.748 


114.4 


216.8 


150.5 


S19 


323.0 


36.140 


0.863 


74.6 


343.7 


157.0 


S21 


35.0 


8.203 


0.784 


55.3 


254.8 


178.0 


S24 


444.0 


44.677 


0.941 


106.6 


3.0 


290.1 


S27 


95.9 


16.072 


0.963 


93.3 


192.1 


312.9 


S29 


85.1 


14.839 


0.936 


127.0 


157.3 


346.3 


S31 


64.2 


12.303 


0.699 


115.1 


138.9 


343.4 


S33 


95.0 


15.967 


0.731 


44.5 


82.7 


326.9 



TABLE 2 

Parameters of the simulations 



Run 


Run 




MiMBH 


a 


P 


e = 0.0 


e = 0.7 




(Mq) 


(mpc) 


(yr) 


1-12 


241-252 


1.0 X 10^ 


400 


0.3 


0.243 


13-24 


253-264 


1.0 X 10^ 


400 


1 


1.481 


25-36 


265-276 


1.0 X 10"^ 


400 


3 


7.696 


37-48 


277-288 


1.0 X 10"^ 


400 


10 


46.837 


49-60 


289-300 


1.0 X 10^ 


400 


30 


243.373 


61-72 


301-312 


2.5 X 10^ 


1000 


0.3 


0.243 


73-84 


313-324 


2.5 X 10-^ 


1000 


1 


1.481 


85-96 


325-336 


2.5 X 10^ 


1000 


3 


7.695 


97-108 


337-348 


2.5 X 10^ 


1000 


10 


46.834 


109-120 


349-360 


2.5 X 10"^ 


1000 


30 


243.355 


121-132 


361-372 


5.0 X 10"^ 


2000 


0.3 


0.243 


133-144 


373-384 


5.0 X 10"^ 


2000 


1 


1.481 


145-156 


385-396 


5.0 X 10"^ 


2000 


3 


7.695 


157-168 


397-408 


5.0 X 10"^ 


2000 


10 


46.828 


169-180 


409-420 


5.0 X 10"^ 


2000 


30 


243.324 


181-192 


421-432 


1.0 X 10-3 


4000 


0.3 


0.243 


193-204 


433-444 


1.0 X 10-3 


4000 


1 


1.480 


205-216 


445-456 


1.0 X 10-3 


4000 


3 


7.692 


217-228 


457-468 


1.0 X 10-3 


4000 


10 


46.816 


229-240 


469-480 


1.0 X 10-3 


4000 


30 


243.263 



The orientations of the massive binary were not selected 
randomly but rather were distributed uniformly on a grid. The 
surface of a sphere can be uniquely tesselated by means of 12 
regular pentagons, the centers of which form a regular do- 
decahedron inscribed in the sphere. We determined the cen- 
ters of the pentagons and transformed the coordinates of the 
points into inclination / and longitude of ascending node fl 
of the IMBH/SMBH orbit. We then rotated the points by an 
arbitrary angle to avoid chance superpositions with the S-star 
orbital inclinations. The resulting distribution of angular mo- 
mentum vectors on the plane of the sky is shown in Figure [TJ 
the values of the 12 angles are given in Table[3] 

We followed the evolution of each 21 -body system for a 
time corresponding to 10 Myr using the AR-CHAIN code 
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Fig. 1. — Directions of the angular momentum vector of the 19 S-stars 
(stars) and of the 12 possible orientations of the IMBH/SMBH orbit (circles). 
In this all-sky map, the vertical dimension refers to the inclination i of the 
orbit and the horizontal dimension to the longitude of the ascending node 
f2, as defined in Eisenhauer et al. (2005) and Gillessen et al. (2009). The 
Galactic plane Ues at / = 0° , H = -3 1 .4° . 



TABLE 3 

Grid of IMBH/SMBH orbital inclinations 



number 


(deg) 


Q. (deg) 


1 


90.0 


31.49 


2 


129.21 


86.17 


3 


50.79 


266.17 


4 


90.0 


211.49 


5 


27.97 


14.14 


6 


81.95 


328.25 


7 


98.05 


148.23 


8 


152.03 


194.14 


9 


66.07 


92.05 


10 


37.09 


169.12 


11 


142.91 


349.12 


12 


113.93 


272.05 




1000 



Fig. 2. — Relevant timescales for the S-stars and the IMBH as a function 
of distance from the SMBH: orbital period (dashed line); star symbols 
are the S-stars, Tprec relativistic precession timescale (solid lines) for eccen- 
tricities e = and e = 0.999, Tqw gravitational wave emission timescale for 
MiMBH = 4OO-4OOOM0 and e = 0-0.7 (dashed region). The vertical dotted 
lines represent the adopted values for the black hole binary initial separation. 
The filled grey region indicates the estimated ages of the S-stars. The vertical 
dotted-dashed line represents the outer edge of the disk of massive stars. 



(accounted for by the PNl and PN2 terms) 
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A^SMBH 
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is generally much shorter than the integration time (see 
also Fig.|2]l. The energy loss due to emission of gravita- 
tional waves (accounted for by the PN2.5 term) occurs on a 
timescale 



(iMikkola & Merrit3l2008h . a recent implementation of the al- 
gorithmic regularization method that is able to reproduce the 
motion of tight binaries for long periods of time with ex- 
tremely high precision. The code com bines the use of the 
chain structure, introduced originally by Mikkola & Aarsedi] 
(U993), with a new time-transformation to avoid singulari- 
ties and achieve high precision for arbitrary mass ratios. In- 
tegrations were carried out on two computer clusters at the 
Center for Computational Relativity and Gravitation at the 
Rochester Institute of Technology: gravitySimulator 
jHarfst et alj |2007|) and newHorizons^, and on one com- 
puter cluster hosted by the Research Computing group at the 
same institution. 

The AR-CHAIN code also includes post-Newtonian correc- 
tions to the accelerations up to order PN2.5. This is necessary 
due to the fact that the general relativistic precession timescale 

^ http://ccrg.rit.edu/facilities 
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where 



and 



F{e)= 



-7/2 



1 73 2 37 4 
1-1- — + — 



24 



A/SMBH * IMBH 



IMBH ■ 



10 ' M, 



Mi, 



(2) 



(3) 



(4) 



MsMBH + MiMBH 

This is typically ^ lOMyr (see Fig.|2]i, implying that the orbit 
of the IMBH does not evolve appreciably over the integration 
time. The only exception is the set of runs with q = 10"^, 
a = 0.3 mpc and e = 0.7, for which 7gw ~ 3 Myr. 
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Some relevant timescales for the S-stars and the IMBH 
are shown in Figure|2] One timescale that is not illustrated 
there is that associated with dynamical friction of the IMBH 
against the background stars. Assumin g a circular orbit for 
the IMBH, the frictional drag is given bv dBinnev & Tremaind 
[l982,Eq.7-18) 



47rlnAGV(r)M 



MBH 



2X 

erf(X)-^e- 
a/tt 



(5) 



where p(r) is the background stellar density, InA the Coulomb 
logarithm andX = v^. / \/2cr is the ratio of the IMBH velocity to 
the one-dimensional velocity dispersion a. The rate of change 
of the angular momentum L = rvc = VGMsmbh r can be writ- 
ten as dL/dt = F r / M\MBn- Assuming X =\, this implies 



dt 



M: 



3/2 



P(r) 



(6) 



SMBH 



The dynamical friction timescale is then 



Tdf = 



1 dr 
r dt 



1 



3/2 
SMBH 



-3/2 



3.47rlnAVGMiMBH p{r) 



2x 10^ yr 



A^SMBH 

4x 10''Mr 



3/2 



lO^M 



Ml 



IMBH 



/lO^Mopc" 

V P 



having taken InA = 6.6 ("S pinnato et al.ll2003h . 

Since our A^-body models include only a fraction of the stel- 
lar mass in this region, the IMBH experiences essentially no 
dynamical friction. This may, or may not, be an accurate 
representation of what would actually happen at the Galac- 
tic center. Recent observations (??) raise questions about 
the existence of a density cusp in the observed stellar pop- 
ulation within the central parsec: Number counts of the late- 
type (old) stars suggest a density that drops inside of ~ 0.5 
pc. Proper-motion data (?) are also consistent with a mass 
density that falls toward SgrA*. Even if a density cusp were 
present initially, an inspiralling IMBH tends to displace the 
background stars and create a low-density core on a scale of 
~ 0.1 pc (Baumgardtet al. 2006). The absence of dynami- 
cal friction in our simulations may therefore be a reasonable 
approximation. 

By ignoring the possibility of a density cusp, we also ex- 
clude the other dynamical effects that would accompany a 
high density of background stars, including cusp-induced pre- 
cession of the orbits; damping of Kozai oscillations; and ejec- 
tion of S-stars by close interactions with stellar remnants. 

3. RESULTS: SHORT-TERM EVOLUTION 

One motivation for basing our A^-body models on the actual 
S-star orbits is that we are able to make quantitative state- 
ments about how these orbits would evolve in the near future. 

While a close encounter of the IMBH with a star can pro- 
duce a sudden large change in the star's velocity, such encoun- 
ters are relatively rare, and the evolution of the star's trajec- 
tory is better described in terms of the change with time of its 
orbital elements. Figures [3] and |4] summarize the evolution of 
the orbit of star SI over the next lOOOyr in each of the 120 
integrations corresponding to ^ = 10"^. We have plotted the 
changes in a, e and Q with respect to their values in the year 
2008; here B is the angle in radians between the star's orbital 
angular momentum vector at the two times. 



The IMBH lies completely inside the orbit of S 1 in the top 
three rows of Figure[3]and|4] In these runs, the star's orbital el- 
ements (which are computed with respect to the instantaneous 
position and velocity of the SMBH) are found to oscillate due 
the movement of the SMBH in response to the IMBH. Oscil- 
lation on both the (short) timescale of the IMBH/SMBH orbit, 
and the (longer) timescale of the star's orbit, can be seen. The 
amplitude of these oscillations is greatly reduced when the or- 
bital elements are computed with respect to the center of mass 
of the IMBH/SMBH binary; there is almost no true evolution 
of the star's orbit on these time scales. 

When the IMBH's orbit is larger than that of the star, oscil- 
latory changes in the star's orbit still take place on the (now 
shorter) timescale of its own orbit, but in addition there are 
stepwise changes that occur roughly once per orbit of the 
IMBH, when it comes closest to the star. The latter changes 
are similar in magnitude for the circular and eccentric bina- 
ries, but in the latter case, changes in Si's orbit are less de- 
pendent on the binary's orientation. 

We chose to place the IMBH at apoastron at the start of the 
integrations, and this choice determined the time of first and 
su^^quent "close encounters" with each star. This is clear 
in the lower panels of Figs. [3]and|4] where the IMBH orbit 
Jiajs/a longer period than that of S 1 . Had we varied the ar- 
— gument of periastron for the IMBH, the jumps in Si's orbital 
pc Elements would have taken place at different times with re- 
sp^g) to 2008. These jumps can be significant, e.g. as large 
as ^ 30 arc minutes in the case of 6. However based on the 
figures, the a priori chance of such a jump occurring during 
any short (^ yea rs-long) interval of ti me is small. 

According to Gillessen et al.l (120091) . the errors in the obser- 
vationally determined values of a, e and 6 for SI are ^ 5%, 
~ 0.03 and ^ 0.01 radians respectively. Figures [3] and |4] 
suggest that short-term few years) changes in the semi- 
major axis might rise above the measurement uncertainties for 
IMBH/SMBH binaries with ^ w 10"^ and a < 1 mpc (e = 0) 
and a < 3 mpc (e = 0.7). Changes in e and might barely be 
detectable on the same timescale for a < 0.3 mpc. 

4. RESULTS: LONG-TERM EVOLUTION 

The perturbations exerted by the IMBH induce changes in 
the orbital elements of the stars. In particular, the presence of 
an IMBH has the following effects: (i) randomization of the 
inclinations of the stars, (ii) ejection of stars from the region, 
(iii) scattering of stars onto plunging orbits that result in tidal 
disruption in the SMBH's tidal field (iv) secular effects like 
Kozai cycles. 

4.1. Randomization of the orbital planes 

The randomizatio n of the orbital planes has been shown by 
iMerritt et al] (120091) to be efficient for IMBH masses larger 
than ~ I5OOM0 and eccentricities larger than ~ 0.5. In 
these cases, a thin corotating disk of stars was converted into 
an isotropic distribution in just ~ 1 Myr. In this work, we 
consider a population of stars that are already approximately 
isotropic at the beginning of the simulations. Nonetheless, the 
IMBH acts to randomize the orbital planes of the stars so that 
the angular momentum vectors can point to very different di- 
rections at the end of the 10 Myr evolution, as can be seen in 
Figure|5]in the case of star S9. 

4.2. Ejections 

Perturbations from the IMBH can also result in ejections 
of stars from the region. While ejections can be produced by 
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Fig. 3. — Changes over the initial 1000 years (i.e. from 2008 to 3008) of the orbital elements of star SI, in runs with q= IQ ,a circular IMBH orbit, and five 
different values of aimbh, from 0.3 mpc (top) to 30mpc (bottom). Different colors denote the 12 different orientations of the IMBH/SMBH orbital plane. 



close encounters with the IMBH, the majority of the events 
are due to slow but steady increase in a star's eccentricity, 
until eventually the star becomes unbound to the SMBH. We 
expect the probability of ejection to depend on both the stellar 
and the IMBH orbital parameters. Figure|6] shows the average 
number of escapers produced in any simulation with given 
mass ratio q = Mimbh/A^smbh and semi-major axis a of the 
binary black hole. We find that most of the ejections occur for 
large mass ratios and for separations comparable to the typi- 
cal distances of the S-stars. These parameters maximize the 
strength of the perturbations, which tend to alter the eccen- 
tricity of the star. A circular and an eccentric binary seem to 
produce approximately the same number of escapers. 

Stars with initially large eccentricity are more susceptible 
to perturbations and therefore more likely to be ejected, as 
shown in Figure]?] In particular, with e = 0.96 star S14 is the 
one which experiences the most ejections. 

Extrapolating the results obtained for the S-stars sample to 
the real Galactic center population (assuming a total of 10"* 
stars in the same region) we obtain an ejection rate 7?. ~ 1 .8 x 

10-4yi-"'. 

4.2.1. Sun'ival time 



study the "survival" probability for the S-stars, for a given set 
of orbital parameters. We estimate the mean time until ejec- 
tion Tej for individual stars for a given set of binary parame- 
ters. Since not all stars are ejected before the end of the sim- 
ulation, an average of the ejection time over the cases where 
ejection occurs would yield a lower limit to Tej. The prob- 
lem of estimating the lifetime of a population given a discrete 
set of events within a finite time is called "survival analysis 
under censoring", where censoring refers to the fact that at 
the end of the simulation only a fraction of the stars have a 
measured time of ejection. We therefore apply a maximum 
likelihood statistical analysis to our sample of censored data 
and compute the mean ejection time for each star and for a 
given set of binary parameters by averaging over the 12 runs 
which correspond to the different binary orientations in the 
sky. We define A^^, the number of ejection events experienced 
by a star by the end of the simulations at time T and Nc the 
number of non ejections. If f, , ; = Q...Ne are the ejection times 
for the events, the mean ejection time for a given star is 



-T 



, ti + —T 

1=0 



(9) 



Since stars can become unbound from the SMBH and es- 
cape the region monitored by observations, it is interesting to 



The results are shown in Figure[8] for all stars and provide 
a quantitative estimate of the probability of ejection for in- 
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Fig. 4.— Like Fig.|3] but for runs with the eccentric IMBH/SMBH binary. 
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Fig. 5. — Orientation of the orbital plane of star S9 as currently observed 
(star) and at the end of the simulations (empty symbols if bound star; full 
symbols if unbound star). Different symbols are for different binary mass 
ratios: ^ = 1.0 X 10^ (circles), § = 2.5 X 10"* (squares), ^ = 5.0 X 10^ (tri- 
angles), 5= 1.0 X 10"^ (diamonds). 



dividual stars taking into account both the orbital properties 
of the stars and of the IMBH. The mean time to ejection 
is ;%10^yr in the case of binaries with x 10"^ and 

3mpc< fl< lOmpc. Similarly to Figure|7] stars with large ec- 
centricities have the shortest Tej . 

4.2.2. Hypervelocity stars 

A fraction of the escapers produced during the evolution 
can be ejected with large enough velocity to become unbound 
to the Galaxy. Such stars would be det ected in the Galac- 
tic halo as hypervelocity stars (HVSs) (iBrown et alJ 12005). 
Assuming a threshold velocity of 700kms~' for escaping the 
Milky Way's potential, we find that about 15% of all the esca- 
pers become hypervelocity stars. The average velocity for the 
full sample of escapers is ^ 350kms~' for the circular binary 
and ^ 400kms~' for the eccentric binary, while for the sam- 
ple of HVSs the mean is ^ 1450kms"' and ^ 1550kms"' 
re spectively. This is con sistent with the theoretical estimate 
bv lYu & Tremaiiig (l2003h 



' 3.2GMsmbhMimbh 

(MsMBH+MiMBH)a 



U \ 1/2 / mpc\ 1/2 / MsMBH 



^2.3 X low (^) (^) 



4x 10''M 



(10) 
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Fig. 6. — Average number of ejected stars as a function of the black hole binary parameters, in the case of an initially circular (left) and eccentric (right) binary. 
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onto plunging orbits, and are therefore rare events in the simu- 
lations. The total number of captures (out of 240 integrations) 
experienced by the stars is 13 and 58 (corresponding to a prob- 
ability of about 5 and 25%) for the circular and eccentric case, 
respectively, indicating that an eccentric binary is more likely 
to lead to a capture event. Large eccentricity stars are more 
likely to be tidally disrupted by the SMBH, as can be seen in 
Figure|9] 

Simulations with q = IQ~^, a = 0.3 mpc and e = 0.7 have 
the shortest time for emission of gravitational waves, namely 
~ 2.7 Myr (see Eq.|2]l. For this choice of parameters, the orbit 
of the IMBH decays due to energy loss, which is correctly 
implemented in the code, until the two black holes coalesce. 
From this point on, the integration is carried on with only one 
black hole. 



Fig. 7. — Fraction of ejections experienced by each S-star as a function 
of the star's initial eccentricity in the case of a circular (stars) and eccentric 
(circles) binary. 

where v = .. — . The distribution of ejection velocities 
for escapers has a peak around ~ 100 kms"' and shows a long 
tail toward high velocities (up to several thousand kms"'). 

4.3. Tidal captures and mergers 

Interactions with the IMBH can also drive stars to be cap- 
tured (either tidally disrupted or accreted) by the massive 
black hole. The Schwarzschild radius for the Milky Way 
black hole is given by 

2GMsMBH 



Rs = 



■3.8 X 10"^ pc 



(12) 



while the tidal radius for a solar type main-sequence star is 

„2,^ \ 1/3 



-3.6x 10-''pc 



/WsMBH 

4o7T06M^ 



1/3 



(13) 



For simplicity we set the capture radius for all stars to R, = 
3Rs. When a star approaches the SMBH within such dis- 
tance, it is merged with the black hole. The mass of the star 
is added to that of the black hole and new positions and ve- 
locities are computed. Mergers occur when stars are scattered 



4.4. Kozai oscillations 

The presence of an IMBH outside the orbits of the S-stars 
results in perturbation s on the inne r stellar orbits. In particu- 
lar, Kozai oscillations (iKozail 19621) can be induced in the stel- 
lar orbits if the SMBH - star - IMBH system can be regarded 
as a hierarchical triple with a large inclination between the in- 
ner and the outer orbit. If the relative inclination is large, the 
time-averaged tidal force exerted on the inner binary produces 
variations in the eccentricity ei„n, inclination j and argument 
of pericenter u. The variations are such that the quantity 



h = il-e^)cosHj), 



(15) 



also known as Kozai integral, is conserved. This is propor- 
tional to the component of the inner angular momentum per- 
pendicular to the outer orbit. The semi-major axis of the inner 
and outer orbits can have periodic perturbations but no long- 
term evolution, so they can be assumed to be roughly con- 
stant. While the evolution of the eccentricity and inclination 
is always oscillatory, the argument of pericenter may either 
oscillate (libration) or circulate {circulation). In either case 
the periodic oscillations in ein,, and j are often called Kozai 
cycles. Remarkably, the variations in /', einn and u do not de- 
pend on the masses of the stars or the dimensions of the orbits. 
These affect only the timescale of the effect and not its size. 

Th e characteristic time-scale for Kozai oscillatio ns is given 
by (iHeggie & Hutjl200l iKinoshita & Nakaill2007h 



4/C 



^out^SMBH + '«* 



3^/611 P'mn /WiMBH 



(1- 



(16) 
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Fig. 8. — Mean time until ejection for all the stars as a function of the binary orbital pai'ameters. 
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Fig. 9. — Total number of mergers experienced by the individual stars as a function of their orbital parameters. Left: initially circular binary. Right: Initially 
eccentric binary. 
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Fig. 10. — Eccentricity evolution of stars S8, S12 and S27 in a run with 
^IMBH = 4OOOM0 , a = 30mpc, e = and with relative inclinations of j = 58° , 
86°, 79°, respectively. Horizontal dashed lines represent the value of the 
maximum eccentricity predicted for the Kozai mechanism while the solid 
vertical lines indicate the Kozai timescale. 

where Pmn and Pout are the period of the inner and outer binary 
and /C is a numerical coefficient which depends only on the 
initial values of j, einn and u. 

The maximum eccentricity e^ax attained by the inner binary 
is also a function of the initial values of j, ei„n and u) only. In 
the case of a circular inner binary, e,nax only depends on the 
initial relative inclination. 

Kozai oscillations can be damped by general relativistic 
precession. 

Three conditions must be satisfied in order for the Kozai 
mechanism to operate in our simulations: (i) the IMBH or- 
bit must lie outside the orbits of the S-stars, so that the three 
body system can be regarded as a hierarchical triple. When- 
ever the inner and outer orbital periods are well separated, 
Eq. [16] accurately predicts the time scale for the oscillations 
in the eccentricity. The amplitude of the oscillations can also 
be predicted from the initial parameters ("Kinoshita & Nakai' 
12007). (ii) The inner and outer orbits must be sufficiently in- 
clined with respect to each other, (iii) The period of the Kozai 
oscillations must be shorter than the age of the system and 
shorter than that of any other mechanism that might wash out 
the oscillations, e.g. relativistic precession. 

We therefore expect to observe Kozai oscillations in the 
simulations with large flout/flinn (like the runs with initial bi- 
nary separation a = 30 mpc) for those stars whose orbits are 
significantly inclined with respect to the SMBH/IMBH plane 
and with a long precession time scale. Also, the timescale for 
oscillations will be shorter in the runs with the largest IMBH 
mass iq = IQT^). 

Figure[TO] shows the evolution of the eccentricity for three 
stars which exhibit Kozai oscillations in a run with Mimbh = 
4OOOM0, a = 30mpc and e = 0. The predicted values of the 
maximum eccentricity and of the Kozai timescale are also 
shown. The agreement is satisfactory in all cases. As ex- 
pected, the semi-major axis of the inner orbit remains constant 



Fig. 11. — Cumulative distribution function for the eccentricities of the 
S-stars after lOMyr of evolution in the case of a circular binary. Different 
lines are for different values of the binary semi-major axis. Each line is an 
average over the 12 different simulations with different IMBH 3D orientation 
with respect to the stars. The four panels refer to the four different values 
of the binary mass ratio. The shaded area represents the initial distribution 
of the stars. The N(e) ^ line con'esponds to an isotropic distiibution of 
eccentricities. 

and h is conserved while the eccentricity and the inclination 
oscillate. 

5. EFFECTS ON THE OBSERVED POPULATION 

The perturbations described in the previous section act to 
modify the orbital elements of the stars. Figure[TT] shows the 
cumulative distribution for the eccentricities of the S-stars that 
remain bound at the end of the simulations. The four pan- 
els refer to the four different values of the binary mass ratio. 
The distributions for Mimbh = 400- IOOOMq are consistent 
with an isotropic distribution of eccentricities P{e) ^ e for all 
values of the binary semi-major axis. For larger black hole 
masses (Mimbh = 2000- 4000 M©), however, only very small 
(a = 0.3 - 1 mpc) or very large {a > 30mpc) semi-major axes 
result in final distributions consistent with observations. Re- 
sults from an eccentric black hole binary appear qualitatively 
similar The figure suggests that the IMBH exerts the largest 
perturbations at separations of order a ~ 3 - 10 mpc, which 
is comparable to the semi-major axes of the S-stars. A close 
inspection of Figure[TT]reveals that the final distributions can 
deviate from the observed one in two ways: either the final 
curve lies above the observed distribution, as is the case for 
a = 1 - 3 mpc, or it lies below the observed distribution, as in 
the case of a = lOmpc. The reason for such distinctively dif- 
ferent behavior is the competition between two effects of the 
IMBH perturbations: the tendency to increase the mean ec- 
centricity of the population and the removal of stars from the 
population once they become unbound. The former process 
tends to move the distribution below the thermal line while 
the latter, by removing the most eccentric stars, tends to push 
the distribution above the thermal line. For semi-major axes 
fl = 1 - 3 mpc, the IMBH lies almost completely within the or- 
bits of the S-stars. Therefore, interactions only occur with the 
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Fig. 12. — Cumulative distribution function for the semi-major axes of the 
S-stars after lOMyr of evolution in the case of a circular binary. Line styles 
are as in Fig. llll 

Stars that have small pericenter distances, leading to ejections. 
This is consistent with the fact that the most eccentric stars 
are the most likely to be ejected (see Fig.|7]i as well as with 
the fact that the largest number of escapers is produced for 
a = 3 mpc, regardless of the binary mass ratio (see Fig.|6]l. For 
a = lOmpc, the IMBH orbit intersects those of all the S-stars, 
thus resulting in strong perturbations on the orbital elements. 
The global increase in the eccentricities wins over the ejection 
of unbound stars in this case. The relative importance of es- 
capers over the global perturbations on the orbital elements is 
also evident if we look at the mean eccentricity of the sample 
for a given binary mass ratio and initial separation. For a = 0.3 
and a = 30 mpc, the value of the mean eccentricity e remains 
close to that of the S-stars, which is ^ 0.75. For a = 1 - 3mpc, 
when ejections dominate, the mean eccentricity is reduced to 
^ 0.6-0.65, while for lOmpc e achieves the largest value of 
-0.8. 

Figure[T2] shows the cumulative distribution of the semi- 
major axes of bound stars at the end of the simulations. We 
limit the range of semi-major axes to 45 mpc, which is the 
maximum value observed in the current S-stars sample. Devi- 
ations from the observed distribution occur preferentially for 
large binary mass ratios and, in these cases, for binary sepa- 
rations in the range 3 - 10 mpc, similarly to what we observe 
for the eccentricity distributions. 

Both FigurefTTI and [T2l show distributions after lOMyr of 
integration. The evolution in the distributions, however, is 
expected to be time-dependent since the perturbations from 
the IMBH and the ejections tend to grow in time (see Fig.O. 
We find that the effect of ejections becomes apparent after 
only 1-2 Myr and after 5 Myr the distributions take the form 
they have at the end of the simulations. 

6. DISCUSSION 

Our results help place limits on the orbital parameters of 
possible IMBHs in the Galactic center Constraints have been 
placed in the past by several authors. One natural constraint 



Fig. 13. — Constraints on the orbital parameters of a hypothetical IMBH in 
the Galactic region. The shaded areas represent regions of parameter space 
that can be excluded based on observational or theoretical arguments. The 
dotted lines mark the distances at which the S-stars are currently observed. 
The dashed line represents the 5 yr orbital period corresponding to discov- 
erable systems. The parameters enclosed in the empty rectangular box are 
required for an efficient randomization of inclinations in the cluster infall 
scenario tMerritt et aUi2009i) . The small rectangular region just below the 
empty box represents the parameter space excluded by this work. 

comes from the requirement that the center of mass of the bi- 
nary coincides with the peak of the stellar distribution within 
the observational uncertainties (Yu & Tremaine 2003.) 



Ml 



IMBH 



MiMBH + MsMBH 



flBH^Smpc, 



(17) 



where cjbh represents the binary semi-major axis. The corre- 
sponding prohibited region in the Mimbh-abh plane is shown 
in the upper right corner of Figure[T3] 

The orbit of a black hole binary shrinks due to emission of 
gravitational waves on a timescale given by Eq.|2] Demand- 
ing a lifetime of at least lOMyr excludes the lowest region in 
F ig.[T3]corresporiding to small separations. 

iGillessen et al.l (120091) derive an upper limit of 0.02Msmbh 
to the mass that can be contained within the orbit of star S2. 
This implies that no more than 8 x 10"* M0 can be hidden be- 
tween — 0.3 mpc and — 9 mpc of the SMBH. 

The most stringent constraints can be derived from the as- 
trometric wobble of the radio source SgrA*, which could 
show indication of binary motion. Very Long Baseline Ar- 
ray observations from 1995 to 2000 show a residual motion 
for SgrA* of — 0.5 mas = 0.02 mpc and an upper limit of 
— 8 km s"' for the peculiar motion perpendicular to the Galac- 
tic plane. Such limits exclude the rightmost triangular region 
in Fig.flJ]) (Hansen & Milosavljevic 2003; Yu & Tremaina 

t03). More recent measurements from IReid & Brunthalen 
)04.) find an upper limit to the motion in Galactic latitude 
of -0.4± 0.9kms"^ Combining this restriction with the 5yr 
limit on observable systems results in the large triangle region 
in Fig.ini 
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The perturbations exerted by an IMBH on the orbits of the 
S-stars allow us to exclude the region of parameter space cor- 
responding to masses 2OOO-4OOOM0 and initial semi-major 
axes ^2 - lOmpc. Such region is represented by the shaded 
box in the figure. 

Interestingly, the IMBH parameters req uired for an efficient 
randomization of inclinations (Merritt et al. 2009) in the clus- 
ter infall scenario (Mimbh ~ 1500Mq for the simulated range 
of separations 10-50mpc; see rectangular box in Fig.[T3]l are 
consistent with all the constraints placed so far Mass ratios 



larger than 10 ^ were not simulated but they are expected to 
operate even faster than the smaller ones. 
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